home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / tos / updates / update27.zoo / pml / diffs
Encoding:
Text File  |  1992-12-28  |  8.2 KB  |  321 lines
  1. ===================================================================
  2. RCS file: /net/acae127/home/bammi/etc/src/master/atari/pml/pmlsrc/Changelog,v
  3. retrieving revision 1.17
  4. diff -c -r1.17 Changelog
  5. *** 1.17    1992/08/14 15:18:43
  6. --- Changelog    1992/12/28 08:11:47
  7. ***************
  8. *** 293,295 ****
  9. --- 293,300 ----
  10.       add target for ident.o
  11.   
  12.   ---------------------------- Patchlevel 19 ------------------------------
  13. + 
  14. + sqrt.c:: Michal Jaegermann
  15. +     major hacks. see the comments at the top.
  16. + 
  17. + ---------------------------- Patchlevel 20 ------------------------------
  18. ===================================================================
  19. RCS file: /net/acae127/home/bammi/etc/src/master/atari/pml/pmlsrc/PatchLev.h,v
  20. retrieving revision 1.16
  21. diff -c -r1.16 PatchLev.h
  22. *** 1.16    1992/08/14 15:20:26
  23. --- PatchLev.h    1992/12/28 08:11:48
  24. ***************
  25. *** 1,4 ****
  26. ! #define    PatchLevel "19"
  27.   
  28.   /*
  29.    *
  30. --- 1,4 ----
  31. ! #define    PatchLevel "20"
  32.   
  33.   /*
  34.    *
  35. ===================================================================
  36. RCS file: /net/acae127/home/bammi/etc/src/master/atari/pml/pmlsrc/math.h,v
  37. retrieving revision 1.14
  38. diff -c -r1.14 math.h
  39. *** 1.14    1992/08/14 15:18:43
  40. --- math.h    1992/12/28 08:11:51
  41. ***************
  42. *** 46,51 ****
  43. --- 46,57 ----
  44.   extern "C" {
  45.   #endif
  46.   
  47. + #ifdef __TURBOC__
  48. + 
  49. + #include <tcmath.h>
  50. + 
  51. + #else
  52. + 
  53.   #ifndef __STRICT_ANSI__
  54.   /*
  55.    *    Create the type "COMPLEX".  This is an obvious extension that I
  56. ***************
  57. *** 76,94 ****
  58.       double        retval; /* val to return */
  59.   };
  60.   
  61. ! #define M_LN2        0.69314718055994530942
  62. ! #define M_PI        3.14159265358979323846
  63. ! #define M_SQRT2        1.41421356237309504880
  64. ! #define M_E        2.7182818284590452354
  65. ! #define M_LOG2E        1.4426950408889634074
  66. ! #define M_LOG10E    0.43429448190325182765
  67. ! #define M_LN10        2.30258509299404568402
  68. ! #define M_PI_2        1.57079632679489661923
  69. ! #define M_PI_4        0.78539816339744830962
  70. ! #define M_1_PI        0.31830988618379067154
  71. ! #define M_2_PI        0.63661977236758134308
  72. ! #define M_2_SQRTPI    1.12837916709551257390
  73. ! #define M_SQRT1_2    0.70710678118654752440
  74.   
  75.   #endif /* __STRICT_ANSI__ */
  76.   
  77. --- 82,100 ----
  78.       double        retval; /* val to return */
  79.   };
  80.   
  81. ! #define M_LN2                0.69314718055994530942
  82. ! #define M_PI         3.14159265358979323846
  83. ! #define M_SQRT2              1.41421356237309504880
  84. ! #define M_E          2.7182818284590452354
  85. ! #define M_LOG2E              1.4426950408889634074
  86. ! #define M_LOG10E     0.43429448190325182765
  87. ! #define M_LN10               2.30258509299404568402
  88. ! #define M_PI_2               1.57079632679489661923
  89. ! #define M_PI_4               0.78539816339744830962
  90. ! #define M_1_PI               0.31830988618379067154
  91. ! #define M_2_PI               0.63661977236758134308
  92. ! #define M_2_SQRTPI   1.12837916709551257390
  93. ! #define M_SQRT1_2    0.70710678118654752440
  94.   
  95.   #endif /* __STRICT_ANSI__ */
  96.   
  97. ***************
  98. *** 183,188 ****
  99. --- 189,196 ----
  100.   __EXTERN double ldexp    __PROTO((double, int));
  101.   __EXTERN double frexp    __PROTO((double, int *));
  102.   #endif /* !_M68881 */
  103. + 
  104. + #endif /* __TURBOC__ */
  105.   
  106.   #ifdef __cplusplus
  107.   }
  108. ===================================================================
  109. RCS file: /net/acae127/home/bammi/etc/src/master/atari/pml/pmlsrc/sqrt.c,v
  110. retrieving revision 1.8
  111. diff -c -r1.8 sqrt.c
  112. *** 1.8    1992/02/03 20:19:23
  113. --- sqrt.c    1992/12/28 08:11:53
  114. ***************
  115. *** 17,23 ****
  116.    ************************************************************************
  117.    */
  118.   
  119. ! 
  120.   /*
  121.    *  FUNCTION
  122.    *
  123. --- 17,23 ----
  124.    ************************************************************************
  125.    */
  126.   
  127. ! 
  128.   /*
  129.    *  FUNCTION
  130.    *
  131. ***************
  132. *** 116,141 ****
  133.    *    machine I tried.
  134.    *
  135.    */
  136. !  
  137.   /*
  138.    *  MODIFICATIONS
  139.    *
  140. !  *    This routine modified to use polynomial, instead of rational,
  141. !  *    approximation with coefficients
  142. !  *
  143. !  *        P0  0.22906994529e+00
  144. !  *        P1  0.13006690496e+01
  145. !  *        P2 -0.90932104982e+00
  146. !  *        P3  0.50104207633e+00
  147. !  *        P4 -0.12146838249e+00
  148. !  *
  149. !  *    as given by Hart (op. cit.) in SQRT 0132.  This approximation
  150. !  *    gives around 5 digits correct.  Two applications of Heron's
  151. !  *    approximation will give more then practically achievable
  152. !  *    13-15 decimal digits.  Multiplications by powers of 2 are
  153. !  *    replaced by explicit calls to ldexp.
  154.    *
  155. !  *    Michal Jaegermann, 24 October 1990
  156.    */
  157.   
  158.   #if !defined (__M68881__) && !defined (sfp004)
  159. --- 116,146 ----
  160.    *    machine I tried.
  161.    *
  162.    */
  163. ! 
  164.   /*
  165.    *  MODIFICATIONS
  166.    *
  167. !  *    Function sqrt(x) is initially approximated on an
  168. !  *      interval [0.5..2.0] by a quadratic polynomial with
  169. !  *    the following coefficients
  170. !  *
  171. !  *         P0  0.3608580479718948e+00
  172. !  *         P1  0.7477707028388739e+00
  173. !  *         P2 -0.1105464728191369e+00
  174. !  *
  175. !  *    These coefficients were computed with a help of Maple
  176. !  *    symbolic algebra system.  Longer approximation interval
  177. !  *    is used in order to avoid multiplication by SQRT2
  178. !  *    constant in case of odd exponents (this idea comes from
  179. !  *    Olaf Flebbe, flebbe@tat.physik.uni-tuebingen.de).  With
  180. !  *    64 bit IEEE representation three Heron's iterations are
  181. !  *    good enough to satisfy essential requirements of Paranoia
  182. !  *    test suite. An absolute error after three iterations is
  183. !  *    below 2e-19 over the whole range (below 5e-10 after two).
  184. !  *    Multiplications by powers of 2 are performed by explicit
  185. !  *    calls to ldexp.
  186.    *
  187. !  *    Michal Jaegermann, 21 October 1992
  188.    */
  189.   
  190.   #if !defined (__M68881__) && !defined (sfp004)
  191. ***************
  192. *** 144,164 ****
  193.   #include <math.h>
  194.   #include "pml.h"
  195.   
  196. ! #ifdef OLD
  197. ! #define P0 0.594604482            /* Approximation coeff    */
  198. ! #define P1 2.54164041            /* Approximation coeff    */
  199. ! #define Q0 2.13725758            /* Approximation coeff    */
  200. ! #define Q1 1.0                /* Approximation coeff    */
  201. ! 
  202. ! #define ITERATIONS 3            /* Number of iterations    */
  203. ! 
  204. ! #endif
  205. ! 
  206. ! #define P0  0.22906994529e+00        /* Hart SQRT 0132 */
  207. ! #define P1  0.13006690496e+01
  208. ! #define P2 -0.90932104982e+00
  209. ! #define P3  0.50104207633e+00
  210. ! #define P4 -0.12146838249e+00
  211.   
  212.   static char funcname[] = "sqrt";
  213.   
  214. --- 149,157 ----
  215.   #include <math.h>
  216.   #include "pml.h"
  217.   
  218. ! #define P0  0.3608580479718948e+00
  219. ! #define P1  0.7477707028388739e+00
  220. ! #define P2 -0.1105464728191369e+00
  221.   
  222.   static char funcname[] = "sqrt";
  223.   
  224. ***************
  225. *** 165,183 ****
  226.   double sqrt (x)
  227.   double x;
  228.   {
  229. - #ifdef OLD
  230.       int k;
  231. -     register int bugfix;
  232. -     register int kmod2;
  233. -     register int count;
  234. -     int exponent;
  235. -     double y;
  236. - #else
  237. -     int k;
  238. - #endif
  239.       double m;
  240.       double u;
  241.       struct exception xcpt;
  242.   
  243.       if (x == 0.0) {
  244.       xcpt.retval = 0.0;
  245. --- 158,168 ----
  246.   double sqrt (x)
  247.   double x;
  248.   {
  249.       int k;
  250.       double m;
  251.       double u;
  252.       struct exception xcpt;
  253. +     register double (*dfunc)(double, int) = ldexp;
  254.   
  255.       if (x == 0.0) {
  256.       xcpt.retval = 0.0;
  257. ***************
  258. *** 192,222 ****
  259.       }
  260.       } else {
  261.       m = frexp (x, &k);
  262. - #ifdef OLD
  263. -     u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
  264. -     for (count = 0, y = u; count < ITERATIONS; count++) {
  265. -         y = 0.5 * (y + (m / y));
  266. -     }
  267. -     if ((kmod2 = (k % 2)) < 0) {
  268. -         y /= SQRT2;
  269. -     } else if (kmod2 > 0) {
  270. -         y *= SQRT2;
  271. -     }
  272. -     bugfix = 2;
  273. -     xcpt.retval = ldexp (y, k/bugfix);
  274. - #else
  275. -     u = (((P4 * m + P3) * m + P2) * m + P1) * m + P0;
  276. -     u = ldexp((u + (m / u)), -1);    /* Heron's iteration */
  277. -     u += m / u;            /* and a part of the second one */
  278.       if (k & 1) {
  279. !         u *= SQRT2;
  280.       }
  281.       /*
  282.        * here we rely on the fact that -3/2 and (-3 >> 1)
  283. !      * do give different results
  284.        */
  285. !     xcpt.retval = ldexp (u, (k >> 1) - 1);
  286. ! #endif
  287.       }
  288.       return (xcpt.retval);
  289.   }
  290. --- 177,206 ----
  291.       }
  292.       } else {
  293.       m = frexp (x, &k);
  294.       if (k & 1) {
  295. !         m = dfunc(m, 1);
  296. !         /*
  297. !          * multiply by 2 if our exponent was odd;
  298. !          * odd k is now too big by 1 but (k >> 1) will take
  299. !          * care of that
  300. !          */
  301.       }
  302. +     u = (P2 * m + P1) * m + P0;    /* the first guess   */
  303. + 
  304. +     u = dfunc((u + (m / u)), -1);    /* Heron's iteration */
  305. + #ifndef __SOZOBON__
  306. +     /*
  307. +      * with current floating point representation used
  308. +      * by Sozobon C we are already below achievable precision
  309. +      */
  310. +      u = dfunc((u + (m / u)), -1);   /* one more          */
  311. + #endif /* __SOZOBON__ */
  312.       /*
  313.        * here we rely on the fact that -3/2 and (-3 >> 1)
  314. !      * do give different results;
  315. !      * the last iteration and adjust final exponent properly
  316.        */
  317. !     xcpt.retval = dfunc (u + (m / u), (k >> 1) - 1);
  318.       }
  319.       return (xcpt.retval);
  320.   }
  321.